#include "solve.h"
#include "matrix.h"

double* gauss(Matrix a, double* p)          //高斯消元法解线性方程组
{
    for (int j = 1; j < a.col; j++)
    {
        for (int i = j+1; i < a.row + 1; i++)
        {
            a.datap[(i - 1) * a.col + j - 1] /= a.datap[(j - 1) * (a.col + 1)];
        }
        for (int k = j+1; k < a.col + 1; k++)
        {
            for (int l = j+1; l < a.row+1; l++)
            {
                a.datap[(l - 1) * a.col + k - 1] -= a.datap[(l - 1) * a.col + j - 1] * a.datap[(j - 1) * a.col + k - 1];
            }
        }
    }
    idown(a, p);
    up(a, p);
    for (int w = 0; w < a.row; w++)
    {
        std::cout << p[w] <<  std::endl;
    }
}

double* lgauss(Matrix a, double* p)         //列主元高斯消元法解线性方程组
{
    int u[a.col];
    memset(u, 0, a.col);
    for (int j = 1; j < a.col; j++)
    {
        double max = fabs(a.datap[(j-1)*(a.col+1)]);
        int index = j;
        for (int m = j+1; m < a.row+1; m++)
        {
            if (fabs(a.datap[(m-1)*a.col+j-1])>max)
            {
                max = fabs(a.datap[(m - 1) * a.col + j - 1]);             
                index = m;
            }
        }
        for (int n = 1; n < a.col+1; n++)
        {
            std::swap(a.datap[(j - 1) * a.col + n - 1], a.datap[(index - 1) * a.col + n - 1]);
        }
        std::swap(p[j - 1], p[index - 1]);
        u[j-1] = index;
        for (int i = j + 1; i < a.row + 1; i++)
        {
            a.datap[(i - 1) * a.col + j - 1] /= a.datap[(j - 1) * (a.col + 1)];
        }
        for (int k = j+1; k < a.col + 1; k++)
        {
            for (int l = j+1; l < a.row+1; l++)
            {
                a.datap[(l - 1) * a.col + k - 1] -= a.datap[(l - 1) * a.col + j - 1] * a.datap[(j - 1) * a.col + k - 1];
            }
        }
    }
    idown(a, p);
    up(a, p);
    for (int w = 0; w < a.row; w++)
    {
        std::cout << p[w] <<  std::endl;
    }
}

/*int main()
{
    double p[9] = {1, 3, 5, 2, 4, 6, 7, 9, 8};
    Matrix a(3, 3, p);
    //double b[3] = {9, 12, 24};
    double c[3] = {14, 20, 47};
    //gauss(a, b);
    lgauss(a, c);
}*/